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Abstract 


An optimization approach to rocket nozzle design, based on computational fluid 
dynamics (CFD) methodology, is investigated for low-Reynolds-number cases. This 
study is undertaken to determine the benefits of this approach over those of classical design 
processes such as Rao’s method. A CFD-based optimization procedure, using the 
parabolized Navier-Stokes (PNS) equations, is used to design conical and contoured 
axisymmetric nozzles. The advantage of this procedure is that it accounts for viscosity 
during the design process; other processes make an approximated boundary-layer 
correction after an inviscid design is created. Results showed significant improvement in 
the nozzle thrust coefficient over that of the baseline case; however, the unusual nozzle 
design necessitates further investigation of the accuracy of the PNS equations for modeling 
expanding flows with thick laminar boundary layers. 

Introduction 

Two low-Reynolds-number rocket nozzles have been designed by using a CFD 
based optimization procedure. The investigation entailed the following five steps: 

(1) Link a well-tested computational dynamics (CFD) code to an optimization 
computer code. 

(2) Create an objective function for the thrust. 

(3) Obtain previously studied nozzle geometry that were classically designed. 

(4) Optimize those cases with CFD. 

(5) And compare the thrust coefficients of the CFD based designs with those of the 
classically designed cases. 

The study was conducted to determine whether CFD contributed a significant performance 
advantage over classical design methods in the optimization of rocket nozzles with large 
boundary layers. CFD-based optimization has proved highly beneficial in the design of 
hypersonic wind-tunnel nozzles, in which the viscous effects are very high. 1 If CFD-based 
optimization proves beneficial for low-Reynolds-number rocket nozzles, then significant 
advancement will be realized in the design process of resistojets, arcjets, and other nozzles 
that are designed for low chamber pressures. 

Relevance and Background 

With an increasing demand for the reduced cost of space endeavors, researchers are 
forced to investigate new technologies, such as improving fuel economy. Greater fuel 
efficiency can be achieved in satellites and spacecraft with the design of better propulsion 
and positioning systems. A push for improved technology in this area has led to the 
development of rocket nozzles that can operate effectively and reliably in space. 
Resistojets, for example, meet the high efficiency and high reliability requirements of space 
travel. 2 These rockets are not intended for the initial launch function, nor are they able to 
create thrust forces comparable to those of conventional rockets; rather, their purpose is to 
produce small amounts of thrust over long periods of time, to gradually accelerate a 
spacecraft up to a certain speed. To meet these requirements, a resistojet nozzle must be 
able to operate at low chamber pressures and must be as small as possible in size. Low- 
Reynolds-number nozzles differ from conventional high-Reynolds-number rocket nozzles 
because of the large boundary layer that is present at the nozzle exit. Despite their 
unconventionality, these low-Reynolds-number nozzles are still primarily designed using 
conventional means. Furthermore, despite the need for efficiency, relatively little effort has 
been devoted to developing methods for designing low-Reynolds-number thrust nozzles. 



Existing Methods for Rocket Nozzle Design 

In the early years of rocket nozzle design, conical nozzles were designed because 
methodology did not exist for the design of efficient contoured nozzles. Although they 
were necessary at the time, conical nozzles were an inadequate solution. Such nozzles 
resulted in performance losses because they did not produce axially directed exhaust; in 
some cases, these nozzles were excessively heavy because of their length. Short and 
efficient contoured nozzles were finally developed when the truncated perfect nozzle 
approach was pioneered; in this approach, a wind-tunnel nozzle was designed for uniform 
flow and then truncated at a much shorter length. Ahlberg et al . 3 present some results that 
were obtained by using this method; Hoffman 4 discusses some variations of the original 
method. 

Within the same time period, an alternate approach to this problem, based on the 
calculus of variations, was formulated. This method was first investigated by Guderley 
and Hantsch 5 and was improved upon by Rao . 6 According to the procedure outlined by 
Rao, the length, ambient pressure, and flow properties in the immediate vicinity of the 
throat are the governing conditions under which thrust is maximized. This classical 
process assumes that an isentropic flow exists in the nozzle; thus, a variational integral is 
formulated by taking into account an appropriately selected control surface. The solution of 
the variational integral yields certain flow properties on the control surface, and the nozzle 
contour is constructed by the method of characteristics to give the desired flow. 
Historically, Rao’s method has been considered the best, as other classical approaches fail 
to surpass its performance. However, the one drawback to Rao’s method, which is its 
dependence on an inviscid design, leaves room for improvement. 

Classical procedures are based on an approximation of viscous effects; these older 
methods rely on an inviscid design (such as the Rao’s method of design). After an inviscid 
design has been completed, a boundary-layer correction is added to compensate for the 
viscous effects. Because, inviscid and viscous effects could not be calculated 
simultaneously, the effects of such frictional phenomenon were approximated with a 
boundary-layer correction. Although this type of procedure is reliable to a certain extent, 
an inherent amount of error exists in any approximation. For instance, in the design of 
hypersonic wind-tunnel nozzles, Candler and Perkins 7 showed that classical methods break 
down when the boundary layer reaches the same order of magnitude as the local radius. 
Similar results were also found by Kim 8 for low-Reynolds-number rocket nozzles. Thus, 
as the Reynolds number of a nozzle decreases so does the accuracy of the boundary-layer 
correction. In turn, Rao’s design method becomes less effective because it relies on the 
boundary-layer approximation to correct for viscous effects. An alternative to calculating 
the inviscid and viscous flows separately is to numerically solve the Navier-Stokes (NS) 
equations; with these equations, viscous effects can be accurately determined. In the past, 
computational capability was such that the NS equations could not be used in the design of 
contours. More recent advances in computational technology have allowed scientists to 
calculate the NS equations, which previously had to be simplified for computation. Such 
advances have given rise to efficient CFD codes that have eliminated the need to 
approximate viscous effects in aerodynamic design. 

In this study, CFD-based design was predicted to show significant improvement 
over classical design procedures. To test this hypothesis, thrust optimization in resisto-jet 
nozzles was chosen because it is an area likely to benefit significantly from CFD-based 
design. Kim 5 states, “ongoing research in low thrust space propulsion has resulted in 
many high-performance space propulsion rockets such as arcjets, resistojets and magneto 
plasma dynamic thrusters.” Resistojet nozzles have very low Reynolds numbers and are 
not well suited to classical design procedures because such thrusters operate at low 
chamber pressures. The dimensions of such nozzles are small, and the Reynolds number 
is low; consequently, the viscous effects are high. Nevertheless, the classical method of 
design is still widely used in determining nozzle wall contour. A design improvement in 
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such nozzles would be a significant development, and CFD is an appropriate tool for 
realizing this goal. 

CFD-Based Optimization 

The use of CFD codes' in design has usually been to provide a direct analysis of the 
performance of existing or proposed designs; CFD codes have been coupled with design 
procedures for a limited number of applications. The most powerful and general CFD 
design procedures are based on solving an optimization problem. A typical design or 
optimization includes the following steps: 

(1) The design requirements are specified. 

(2) An objective function is constructed, the minimum or maximum of which 
yields the design requirements. 

(3) The set of design parameters or variables is specified. 

(4) An initial value for each of the design parameters is estimated. 

(5) An initial CFD solution is computed by using the estimated design parameters. 

(6) The objective function is computed from the difference between the design 
requirements and the computed solution. 

(7) The sensitivity of the objective function to the design parameters is calculated 
(sensitivity coefficients). 

(8) An optimization problem is solved to generate a new set of design 
variables. 

(9) A new CFD solution is computed and compared with the design requirements. 

(10) If the design requirement is met or a minimum or maximum is reached, then 
the procedure stops, otherwise the process is repeated from step 6 onward. 

In this investigation, the above-mentioned optimization procedure has been applied 
to the design of low-Reynolds-number rocket nozzles. For computational efficiency, the 
nozzle flow field was computed in two steps. The first step was to compute the subsonic 
and transonic flow regions with the NS equations. This step is time iterative, which makes 
it a time consuming but necessary step. Upon convergence of the NS equations, a subset 
of the same equations can be used to more quickly determine the supersonic portion of the 
flow field. The subset of formulas is known as the parabolized Navier-Stokes (PNS) 
equations. These formulas are simplified or parabolized by neglecting both streamwise 
diffusion effects and a portion of the subsonic streamwise pressure gradient. As a result, 
the PNS equations can be integrated by using efficient space-marching techniques. For 
many practical flow fields, the approximations made by the PNS code are valid as long as 
streamwise separation does not occur or as long as the expansion is not too large. For a 
more detailed description of this procedure, see Korte et al. 

With the increased use of supercomputers, both CFD and CFD-based design 
optimization have become significantly faster and more efficient. In the current study these 
new techniques were successfully used to improve the existing design process for low- 
Reynolds-number rocket nozzles. Gains in thrust of greater than 3 percent have been 
obtained for a very simple nozzle wall geometry. This paper discusses the results and the 
methodology used. In this work, the CFD-based thrust optimization process is 
demonstrated by optimizing a low-Reynolds-number thrust nozzle for conditions that have 
been previously investigated by other researchers. 

Numerical Methods 

The CFD solutions for the nozzles are calculated with the code described in Ref. 9. The 
nozzle is divided into two regions for improving computational efficiency: the subsonic- 
transonic section and the supersonic-hypersonic section. In the subsonic-transonic section, 
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the fall NS equations are solved. In the supersonic-hypersonic section, the PNS equations 
are solved. 

Transformation to Computational Coordinates 

The physical domain is transformed into computational space by using the generalized 
transformation 


%(x,y)=& T](x,y)=rj (1) 

where £ represents the streamwise direction and r\ is the crossflow direction. 

The nozzle wall radius or height y w ( x) is defined by using a set of points, the 
coefficients to a spline, or an analytical expression. For a given streamwise location, the 
stretching transformation of Roberts can be used either to cluster points at the wall only or 
at the wall and the centerline. 10 For the extremely thin boundary layers typically found in 
the throat region, neither of these stretching functions is satisfactory. If points are clustered 
only near the wall, then an insufficient number of points is at the centerline and the solution 
deteriorates in that region. Likewise, if points are stretched at both the wall and the 
centerline, then an insufficient number of points is located in the middle of the domain to 
accurately compute the flow field. To overcome this problem, two different stretching 
functions can be blended together by using a fourth-order polynomial in rj . 1 1 
This transformation is given by 

y(j)/yw(x) = 

(i - 7 j 4 >(i - a + 2 a /{i + [(& + 1 m - Df-n) 

_ ( 2 ) 

+rf - 2/5 w /{ 1 + [(ft, + l)/(p w - l)f }) 

n = O'- -i) 'Uw-D (3) 


The coefficient is a positive number greater than 1. The closer that ft is to 1, the more 
clustered the grid becomes. The first line of equation (2) clusters the points at the centerline 

based on the value of , and the second line clusters the points at the wall based on fiw. 

The clustering at the centerline is held fixed at fi c =1 .02, the clustering at the wall is selected 
for each particular case. 

Optimization Setup 

In this section, the assumption is made that the shape of the subsonic-transonic nozzle 
contours of the nozzle has been specified and that an NS solution has been computed. The 
PNS solver used in reference. (9), (with some minor modifications) is coupled to the 
optimization program described below. 

Objective Function 

A nonlinear optimization problem is solved to determine the design parameters by the 
minimization or maximization of an objective function. An objective function Obj is 
dependent on a set of design parameters X: 
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( 4 ) 


Obj(X) = J (pu 2 + p)dA 

o 


where p, u, p, and A are the local density, axial velocity, pressure, and cross-sectional 
area, respectively. For a rocket nozzle, Obj(X) is defined as the exit thrust (into a vacuum) 
which must be maximized to obtain an optimum design. 

Design Parameters 

Like the objective function, the selection of the design parameters depends on the 
problem to be solved. The ideal set would contain the minimum number of elements and 
would be strongly coupled to the objective function. The design parameters are usually 
coefficients used to define a wall boundary or quantities that specify the flow-field 
conditions. For a nozzle design, the flow conditions are usually given, and the wall 
contour must be determined. 

Wall Contour 

Three different functions (i.e. linear, cubic, and sixth-order polynomial equations) are 
used to describe different nozzle wall shapes downstream of the throat. The linear function 
describes a conical nozzle, for which the final slope (or angle) is the design parameter. The 
length is then determined by using the area ratio and the specified angle. For a contoured 
nozzle, defining the wall becomes more complicated. The simplest method in this case 
would be to fit a third-order polynomial to the beginning radius, initial slope, exit slope, 
and length. 


Subsonic-Transonic Section Geometry 

The subsonic-transonic section is specified with a circle: 

x 2 + (y w -R c -r*) 2 = R 2 c ( 5 ) 

where R c and r* are the throat radius of curvature and the throat radius, respectively. 


Linear (Conical Nozzle) Geometry 

The conical nozzle begins at the axial location at which the slope of the subsonic-transonic 
section equals the tangent of the cone angle: 

y w =y 0 + tan ^ (Ax) (6) 

where 6 C is the cone angle. Ax = x - x 0 and x 0 and y 0 are the initial points for the conic 
section. 


Contoured Nozzle Geometry-Cubic Polynomial 

The cubic polynomial begins at the first section used in the PNS calculations and the radius 



y w = y„ + y' 0 (&x) + cAx 2 + dAx 3 
c = 3(y L -y c - y’ 0 L)l(L 2 )-{y' L -y' 0 )l L 
d = Kyi -y:)/L 2 ] + 2(-y L + y o + y' o L) /(L 3 ) 


( 7 ) 


where Ax = x - x 0 , and (x^ y & y 0 ’) is the initial point and slope for the supersonic 
section, (y D y L ’) is the radius and slope at the exit, and L is the length of the contoured 
section. 

The conical and contour nozzle geometry used in this study are shown in Fig. 1. 

Low-Reynolds-Number Thrust Nozzle Design and Optimization 

The CFD-based thrust optimization process was conducted by: 

(1) choosing an existing nozzle design that was previously calculated by other 
researchers; 

(2) varying the grid density of the solution to determine if the specified resolution 
produces accurate results; 

(3) optimizing those cases by using a CFD code coupled to an optimizer, 

(4) and looking for improvements in the thrust coefficient for each case. 

This four-step process was followed to determine the advantages of CFD-based 
optimization over existing methods of design for low-Reynolds-number rocket nozzles. 
Kim 8 investigated 20° and 30° degree conical resistojet nozzles and found that they have 
higher thrust coefficients than contoured resistojet nozzles designed with Rao’s method for 
the same area ratio. For this reason, the optimum angle for a conical resistojet nozzle with 
an area ratio of 82 was determined to provide a target for the CFD-based optimization 
process; thus, if the CFD-based optimization procedure produces a contoured design that 
has a higher thrust coefficient than is obtained from the conical nozzle with the optimum 
angle, the hypothesis is supported and the CFD-based design process is a better design 
process for the studied cases. 

Solution Accuracy and Grid Dependence 

Initial calculations were made with the flow and geometry conditions given by Kim 8 for 
20° and 30° conical nozzles. Viscous flow solutions were computed for hydrogen-gas 
flows through the nozzles. The following assumptions were made: the flow was laminar; 
the back pressure was 0; the wall was adiabatic; the stagnation pressure was 150000 Pa, 
and the stagnation temperature was 1500K. The computed nozzles had a throat radius of 
4.2x1 0‘ 4 m and an exit area ratio of 82. 

Obviously, grid density was an important factor in determining both the precision of a 
solution and die computational time. To minimize the computational time in demonstrating 
the accuracy of the solution, the 20° conical nozzle was run with different normal grid 
densities of 63 and 125 points. The stability requirements of the explicit PNS space- 
marching algorithm force die streamwise grid to have a large number of stations (> 1000), 
which ensures accurate resolution in the streamwise direction. A constant grid-stretching 

factor of p w =1.2 was used in all calculations. Figure 2 compares the streamwise 
momentum profiles at the nozzle exit for the different normal grid densities. Because the 
streamwise momentum is integrated to obtain the nozzle thrust, the profiles are believed to 
be a good measure of the solution accuracy. Very good agreement was obtained for the 
two grid densities; the 125 point solution predicted a slightly lower streamwise momentum 
in the inviscid core. The nozzle thrust coefficients for the 63 and the 125 normal point 
distributions were 1.482 and 1.480, respectively, which is a difference of only 0.14 
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percent. This agreement was sufficient to justify the use of a grid with 63 normal points 
for the nozzle design and analysis. A complete NS calculation was then computed for the 
nozzle to determine the effect of the PNS approximation on the thrust coefficient. A NS 
solution was computed with the same distribution of 63 normal points and 235 axial 
stations. The thrust coefficient for the NS calculation was 1.492, which compares 
favorably with Kim’s 8 NS calculation of 1.498 (603 lb f -sec/lb m vacuum specific impulse) 
with 60 normal grid points. The difference in the thrust coefficients for the NS and PNS 
solutions is less than 0.7 percent. We expect the PNS solutions to follow the same trends 
as the NS solutions, which justifies the use of the PNS equations in the design procedure. 

Low-Reynolds-Number Conical Nozzle Design 

Next, the optimum angle was determined for a low-Reynolds-number conical thrust 
nozzle with a fixed area ratio of 82. The throat geometry was constructed as a circular arc 
(fig. 1) and the PNS calculations began at 1.747x1 O' 4 m downstream of the throat, where 
the wall radius and slope are 4.447x1 O' 4 m and 0.2887, respectively. The subsonic flow 
field of the nozzle was predicted by solving the NS equations to determine the initial 
conditions for the PNS calculations. The PNS code was then modified to construct conical 
geometry as a function of the cone angle (Eq. 6) and to compute the nozzle thrust (Eq. 4). 
The solution had a grid density of 63 normal points, approximately 40,000 streamwise 
stations, and a grid-stretching factor of 1.2. The remaining calculations presented were 
computed on similar grids. 

After the modification was completed, the CFD code was linked to the OptdesX 12 
software system, which is an interactive computer program for computer-aided 
optimization and design. An optimization problem was defined to maximize thrust by 
determining the optimum expansion angle. For each case, the software’s default 
optimization scheme was used (Generalized Reduced Gradient algorithm). The derivative 
of the nozzle thrust with respect to the cone angle was computed with the forward- 
difference selected by OptdesX as the best method with a step size of 0.001. Two cases 
were run using starting angles of 20° and 30°. Because of the simplicity of the objective 
function (i.e. one design parameter), the optimizer required few iterations to achieve the 
optimum expansion angle of 25.96° resulting in a thrust coefficient of 1.488. This result 
was achieved using the different initial conditions for the conical angle. The Mach contours 
for the optimum conical nozzle are shown in Fig. 3. A solution was reached on an Iris 
Indigo after 3 iterations and 25 analysis calls (approximately 3 hr). 

Low-Reynolds-Number Contoured Nozzle Design 

For this problem, the wall contour was defined with a third-order polynomial. Only 
two variables, the length and the exit slope, were varied to determine the optimum 
geometry. A cubic was fitted to the previously determined optimum conical nozzle and was 
used as the starting point for the contour optimization. Then, the initial case was calculated 
and the optimizer varied the parameters in each iteration to generate a new wall contour. 
The CFD code produced a thrust coefficient value for each contour, and the process 
continued until an optimum was reached. The derivatives were computed with a forward 
difference-method with a step size of 0.001. The initial design parameters were specified 
from the optimum conical nozzle. The solution required 23 analysis calls (approximately 3 
hr on an Iris Indigo.) The exit slope converged to the side constraint value of 0.0, and the 
nozzle length (from the throat to the exit) converged to 0.00767 cm. The final thrust 
coefficient was 1.496. The Mach contours for the optimum third-order wall contour are 
given in Fig. 4. 

To determine the effect of the exit angle side constraint on the solution, it was changed 
to -0.364 (-20 degrees) and the problem was restarted. The exit slope converged to the 
side constraint, the length to 0.00779 cm. and the thrust coefficient increase to 1.538. This 
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resulted in an unusual shaped nozzle (Fig. 5). To determine if this was an artifact of 
maintaining the area ratio at a fix value of 82, it was added as a design variable in the 
problem. This solution did not change when the area ratio was included as a design 
variable. The converging nozzle exit (-20 degrees) design is an unexpected result and does 
not follow normal convention. Previously, it had been assumed that the use of the PNS 
solutions in the design procedure would follow the same trends as NS solutions. 
However, the effect of ignoring part of the subsonic streamwise pressure gradient in the 
PNS equations may be more significant in large laminar boundary-layers present in this 
nozzle flow field. Additional investigations comparing the PNS and NS solution for this 
no zzl e design need to be done to determine the validity of this approach and to quantify the 
effect of the approximation made using the PNS solutions for low Reynolds number 
expanding flows. 


Concluding Remarks 

A computational fluid dynamics approach to low-Reynolds-number rocket nozzle 
design has been investigated and compared with previously developed methods. Earlier 
studies have found that conical nozzles produce higher thrust than contoured nozzles 
designed with classical methods. The current study demonstrates a method that is an 
improvement to both the classical method and the conical design based on solutions to the 
parabolized Navier-Stokes (PNS) equations. A CFD code was modified to provide a thrust 
coefficient value for individual cases. The modified analysis software was linked to the 
OptdesX optimization software, and a contoured nozzle design was generated by solving 
an optimization problem. The CFD-based optimization procedure achieved better results 
than those obtained from Rao’s method or from conical nozzles, however, the contoured 
nozzle designs were not as expected. The best results were obtained when the nozzle exit 
angle was allowed to be negative. Further study was suggested to compare the PNS 
results with Navier-Stokes solutions to see if the design parameters for maximum thrust 
coefficient agree and to quantify the accuracy of the PNS solutions for low Reynolds 
number expanding flow fields. 
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(b) Contoured nozzle geometry. 

Fig. 1. Conical and contoured nozzle geometry design parameters. 
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Fig. 2. Axial momentum profile at nozzle exit. 



Fig. 3. Mach number contours for optimum conical nozzle design. 
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